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Abstract. The low temperature properties of double exchange model in triangular 
lattice are investigated via truncated polynomial expansion method (TPEM), which 
reduces the computational complexity and enables parallel computation. We found 
that for the half-filling case a stable 120° spin configuration phase occurs owing to the 
frustration of triangular lattice and is further stabilized by antiferromagnetic (AF) su- 
perexchange interaction , while a transition between a stable ferromagnetic (FM) phase 
and a unique flux phase with small finite-size effect is induced by AF superexchange 
interaction for the quarter-filling case. 
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1 Introduction 

Doped manganite has become one of the most important strongly correlated systems, 
since colossal magnetoresistance (CMR) effect was discovered in 1990s [1}[2|. CMR is re- 
ferred to the resistivity of material change orders of magnitude under external magnetic 
field and it may have a potential application in computer technology or even spintronics. 
There are rich phase diagrams and many ordered phase 0, rising from delicate inter- 
action between electron, spin and orbit degree of freedoms. Further it is found that the 
phase separation [4] (PS) may be crucial to CMR effect. 

Double-exchange model, as a starting point to study manganite, describes the Hund 
interaction between itinerant electron and localized spin of Mn atoms and is expressed 
by 

H DE = -t £ (ClC^+hcj-jH^Clcr^C^-Si, (1.1) 

<ij>,a T-A,f> 

where the nearest-neighbor hopping integral t is adopted as the energy unit, Jh is the 
Hund interaction strength, Cf u (Q A ) creates (annihilates) one electron at site i with spin 
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ol, < ij > stands for the nearest neighbors of lattice site, and cr a c is the Pauli matrix. The 
localized spin S, at site i is assumed as 1 here. In addition, antiferromagnetic (AF) su- 
perexchang is crucial to stabilize AF phase of some underdopped narrow-band mangan- 
ite, and this interaction is expressed by Haf = jAFH<ij>Si-Sj. So the total hamiltonian 
is H = Hjje + Haf- In this model, localized spins is treated as a classic field cp and elec- 
tron degree of freedom can be integrated for any given localized spin configuration. The 
partition function is expressed by Z = Ti c Tr-p(e^^ H ^^l in ^) = Tr c e~ Seff W, and the effec- 
tive action is S e ff(<£) = — J^ v ln(l-\-e~P( £v ~^)+pE((p). Here e v is the i/-th eigenvalue of 
one-electron sector's Hamiltionian matrix H((p), E((p) is the interaction energy between 
spins, /3 is the inverse temperature, ]i is the chemical potential, and n e is the electron 
number operator, respectively. The fluctuation of localized spins is suitable for Monte 
Carlo (MC) simulation, and the partition sum is replaced by stochastic samplings with 
the Boltzmann weight e~ Se "^ / Z. With MC simulation, an intrinsic PS between high 
electron density and low electron density has been reproduced in doped maganites 0. 

Despite the success of reproducing PS, the above method is suffered from finite-size 
effect, since the computational complexity of exact diagonalizing (DIAG) H(cp) scales as 
0(N 4 ), with N being the system size. In order to overcome the above restriction, Fu- 
rukawa and Motome ©[7) proposed Chebyshev polynomials expansion method(PEM) 
and the computational complexity became 0(MN 3 ) at a finite cutoff M. In 2004 they fur- 
ther reduced the computational complexity to O(N) via truncated polynomial expansion 
method (TPEM) |5H9|. The system size HO) is extended to 40 x 40, compared with 10 x 10 
via DIAG [5J. Not only the computational complexity is greatly reduced, but also parallel 
computation is capable under PEM and TPEM, which would increase the computational 
speed very much. Finally, from the viewpoint of physical properties, one-body quantity 
such as energy and electron density, but also dynamical quantity or two-body quantity 
such as conductance are able to be investigated under PEM and TPEM, and there have 
been several examples of application to large scale system |T0Hl3l . 

Up to date, most studies focus on the double-exchange model in one- and two-dimensional 
square lattice. In this paper, we present our results on this model in two dimensional tri- 
angular lattice. One famous layered triangular lattice is the superconductor Nao^Co02- 
I.35H2O discovered in 2003 |14|. Though triangular lattice structure is not found in mag- 
anite yet, it is still meaningful to study the interplay between electron, spin and lattice for 
double exchange model in two dimensional triangular lattice. With the frustration, there 
is 120° spin configuration in low temperature for the half-filling case, and this phase be- 
comes more stable by AF superexchange interaction; but for the quarter-filling case, we 
found a stable long-range ferromagnetic (FM) ordered phase turns into a unique flux 
phase induced by AF superexchange interaction. This paper is arranged as follows: we 
briefly describe PEM and TPEM (more details please refer to Ref. (SHU) and give some 
benchmark in Sec 2, in Sec 3 show results in triangular lattice, and finally we summarize 
in Sec 4. 
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2 The method and benchmark 

The core of Chebyshev polynomial expansion method (PEM) is that the mean value of 
physical quantity is expressed by the sum of the moment multiplied with physical quan- 
tity expansion coefficients, and it converges to exact value guaranteed by rapidly de- 
caying expansion coefficients. It is important that the accuracy is controlled only by the 
cutoff M, the thresholds e p and et r , not by the system size and/ or the chemical potential. 
It is further verified that the error does not accumulate with Monte Carlo sweep. It is 
advantageous that the application to large-scale system is capable, as the computational 
speed is greatly improved due to reduced computational complexity. In following, we 
will give a brief description of PEM and TPEM, and some benchmark. 



2.1 Polynomial Expansion Method 

The Chebyshev polynomials T m (x) are recursively defined by: Tq(x) = 1, 7\(x) = x, and 
T m (x) = 2xT m _i(x) — T m -i{x), with —1 < x < 1. These polynomials are orthonormal in a 
form that 

r-1 dx 

-;Tm(x)T m ,(x)=tx m S mm ,, 



i-i VT^x 1 

{ 1 , m = 0, 

where a m = < i ' ' The density of states (DOS) can be expressed as D(e) = 

2;rv / 1 lje ; Em=oFM^( e )' with the moment being }i m = j\ l T m (e)D{e)de. Then the mean 
value of any physical operator / can be expressed by 

</>= f deD{e)f(e)=Y J }imfm, (2.1) 

J— l m 

with expansion coefficients f m = ^- J_ t -j== x2 f(x)T m (x), which decay exponentially for 

m»l and ensure the accuracy at finite cutoff M \7\. Here f(x) may be the effective 
action function S(x) = — ln[l+e^*~^] or electron number function n{x) = l/[l+e^ x ^f l '] r 
and f m is easily calculated. 

The moment }i m = Y^,=\^m{^v) = Tr(T m (H)) can be evaluated under any set of or- 
thonormal basis e(k), 

Ndim N dim 

Hm=E <e(k)\T m {H)\e{k)>=£ Vm (k), (2.2) 
k=l fc=l 

Here we define a vector v(k,m)=T m (H) \e(k) > and a partial moment ]i m (k) = <e(k)\v(k,m)> . 
From the definition, it is easy to calculate the vector v(k,m) recursively and the element i 
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in the vector v(k,m) is expressed as 

v i (k,0) = e i (k), 

v i {k,l)=Y J H ij Vj{k,Q), 

i 

Vi(k,m)=2j2HijVj(k,m-l)-Vi(k,m-2). (2.3) 

7 

2.2 Truncation of Matrix-vector Product and Trace Operator 

For a sparse hamiltonian matrix Hy with 0(N) nonzero elements, the multiplication can 
be confined to nonzero elements that scale as 0(M D ) where D is the system dimension. 
The computational complexity of M-th moment ^(fc) is 0(M D+1 ) and further reduced 
to 0(M D/2+1 ) if elements less than the threshold e v are neglected. A subspace N £p (k,m) 
is defined by N £p {k,m) = {f\ =Q {i} ,\vi (k,m')\ > e p . The total error of the Boltzmann weight 
isO(M D+1 e p ). 

Based on the importance sampling method, during each MC sweep a new configu- 
ration is accepted as the ratio r = exp~ ASeff is larger than a random number uniformly 
distributed between and 1. The update of the effective action AS e g = Ymk \}i*m N W ~~ 
y°~{Jk)\S m is expressed by 

Ndim 

AS eff = £S m £ Ap m (k), (2.4) 

m k=l 

and the update of the moment is 

Ap m (k) =<e(k) \v new (k,m) >-< e(k)\v old (k,m) > . (2.5) 

Due to local update of the classical field, only a few elements of the Hamiltonian ma- 
trix H are modulated, there are only a limited amount of v(k,m) to be updated, and a 
new threshold e tr can be taken to reduce the computational complexity. The total com- 
putational complexity of one MC sweep is 0(M D+1 N). As Eq. (I2.2D is shown, each basis 
is inter-independent and the calculation of the moment }i m can be parallelized. 

2.3 Benchmark for double exchange model 

Here we present some benchmark, such as the accuracy of the physical quantity, the 
number of nonzero elements, and the sweep time. One basic factor of the algorithm is 
the accuracy. As we mention above, the error can be controlled systematically by M, e p 
and e tr , and the result is often reliable at finite M. In Fig. Q] we show the error of the 
effective action S e ff and the electron occupation n e under a random spin configuration in 
two dimensional triangle lattice. As M increases, the error is smaller and smaller, and the 
result eventually converges to the exact value. From Fig. Qla-b), both the effective action 
and the electron density are very accurate for 20 < M < 40, though M ~ 40 is required at 



5 



lower temperature due to slowly decaying f m . Fortunately the error does not depend on 
other physical parameters, such as the chemical potential ji (Fig. Qlc-d)) and the system 
size (Fig. [H e ~f))- Moreover, the error of the effective action does not accumulate during 
Monte Carlo update process (not shown here). These facts show that PEM combined 
with MC simulation is reliable to study large-scale system. 




Figure 1: The error of S e ff (a, c, e) and n e (b, d, f) vs the cutoff M for a given spin configuration in triangular 
lattice. /h = 8 and /af = 0. Other parameters are (a) and (b) N = 6x6, }l = —S\ (c) and (d) N = 6x6, /3 = 50 ; 
and (e) and (f) /5 = 50, }i = -8. 



The other important factor of the algorithm is the computation speed controlled by 
computational complexity. For TPEM, the complexity is usually associated with N, the 
system dimension D, M, e p and e tr . In Fig. I2|a), we show how the number of elements in 
the subspace N £ (N e ) change with N. While not surprisingly N e decreases with higher 
e p , it is interesting that N e saturates no matter how large N is. That is why the complexity 
is linear with N for one MC sweep. In Fig. |2fb-c) the asymptotic behavior of N e changes 
from 0(M D ) to 0(M D ^ 2 ) once the truncation e p is taken. Now we compare the com- 
putational speed. Fig. [3] illustrates the CPU time cost by one MC sweep via DIAG and 
TPEM, respectively. In one dimension TPEM is prior as N > 64, and it is 1000 times faster 
than DIAG at N > 512. In two dimension, TPEM is advantageous as the system is larger 
than 14 x 14. The sweep time is the same in square lattice and triangular lattice by DIAG 
, but it is shorter in square lattice than in triangular lattice by TPEM, as the hamiltonian 
is sparser owning to electron hopping between less nearest-neighbors in square lattice. 
Due to the dimensionality, the scaling behavior of time with respect to the system size is 
the same in triangular lattice as in square lattice. Once parallel computation is realized, 
TPEM will be more efficient. 
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Figure 2: N e elements in the subspace N e vs (a) the system size N in one dimension, and vs (b, c) the cutoff 
M in one and two dimension. £ = 75, / H = 8, f* = -8, }af = 0- (a) M = 16, (b) N = 40, and (c) N = 40x40. 

3 Double Exchange Model in Triangular Lattice 

In triangular lattice with the lattice constant a, two basic lattice vectors are flj = (l,0)a 
and fl 2 = (i/^)fl. The reciprocal lattices are h = -4^(^,-1) and fc 2 = -^(0,1). For L x L 

triangular lattice, the momentum q = j-b\ + ^b 2 is shortened as ((71,(72)/ with = j^,qi = \, 
and m, n as integers from to L. Here we mainly focus on the spin structure factor S(q) 
defined by 

S(q) = -^E< S rSj>^ 7,i > (3-1) 

y 

where < S ; -Sy > is the mean spin-spin correlation, r« is the displacement between site i 
and j, and N = LxL. 

In order to understand the interplay between electron, spin and lattice, we study the 
spin-spin correlation for different filling case. First we will check the accuracy of the 
mean spin-spin correlation. The maximal monte carlo step is 10000, and the physical 
quantity is evaluated every 20 steps after first 2000 warmup steps. The parameters are 
M = 30, e p = 10~ 5 and et r = 10~ 3 . Due to periodic boundary condition (PBC), we only 
show < Si • Sj > with 1 < j < N in FigSJ For half-filling case (Fig. St a, d)) and low-filling 
case (Fig. H^c, f)) regions, the accuracy are enough. But for quarter-filling case (Fig. Hftb)), 
the accuracy is not good, and higher M is required. For the spin flux phase, M = 30 is 
enough for the accuracy as Fig. 4(e) shows, and the accuracy was verified in 6 x 6 to 
12 x 12 triangular lattice. Therefore it is not necessary to use large M as the basic physical 
phenomenon still remains. 

3.1 Half-filling case 

For half-filling case one lattice has one electron and electron hops between adjacent lat- 
tices will reduce the total energy. According to Pauli exclusion rule two electrons on one 
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Figure 3: Comparison of one MC sweep CPU time (sec) via TPEM and DIAG in (a) one and (b) two dimension. 

M = 30, e„ r = lCT 5 and e„. = l(r 3 . 



lattice site should have opposite spins, hence the nearest-neighbor localized spins tend to 
be antiparallel, but there exists frustration effect induced by triangular lattice. Therefore 
120° spin configuration occurs, where three localized spins of each triangle are coplaner 
with angle between any two being 120°. Using the method in Ref. fl6f, the coplane is 
verified by the fact that (Sj x SA -S^ statistically equals to zero, where site i, j and k are 
located in one triangle; and the angle between any two being 120° is verified by the spin- 
spin correlation of adjacent sites being —0.48, very close to the limit —0.5. Regarding 
to S(q), two peaks are located at (|,|) and (§,§), respectively, and S(q)/N = 0.408 (Fig. 
|5](a)) . This phase is stabilized by AF superexchange interaction with S(q) /N = 0.464(Fig. 
Etd)). 



3.2 Low- and Mediate-filling case 

For low- and mediate-filling case, there is less than one electron on each lattice site 
and adjacent spins tend to be parallel. Hence the system is in FM phase and the peak 
S(0,0)/N = 0.821 is close to the maximum value 1 (Fig. IS^c)). But in the presence of 
superexchange interaction, the system is paramagnetic (Fig. (5]^f )) because of the compe- 
tition between electron-mediated FM correlation and AF superexchange interaction. The 
spin structure factor is almost flat at all vector q and has no peak at all. But for near zero 
filling case, 120° spin configuration will be induced by AF interaction with the interplay 
between spins and lattice. 
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Figure 4: The spin-spin correlation <S\-S;> of 6x6 triangular lattice. /3 = 75 and Jh = 8. £p = 10 and 
ef( . = lCT 3 . (a) J AF =0, ^ = -6 and <«>=0.9272, (b) } AF = 0, ^ = -8 and <n>=0.4512, (c) J AF = 0, ]4 = -10 
and <n>=0.1944, (d) /af=0.1, ?/ = -6and <n>=0.9146, (e) J AF =0.1, f( = -8and <n>=0.5, (f) J AF =0.1, 
fi = -W and <n>=0.2101. 

3.3 Quarter-filling case 

For quarter-filling case in one- and two-dimensional square lattice with PBC, the FM 
phase is not stable, that is positive spin-spin correlation at short distance and negative 
at long distance. This phenomenon has been observed long before (Ref. |5| and the Ref- 
erence inside). Open boundary condition (OBC) or other kind of boundary condition 
is chosen to stabilize FM phase [5J. In two-dimensional triangular lattice with PBC, the 
effect of lattice on FM phase is significant. Each site has 6 nearest-neighbors, and their 
spins tend to be parallel. It is found that electron-mediated FM phase is stablized even 
at long distance and S(0,0)/N = 0.899(Fig. Efb)). So the FM phase is strengthened by 
triangular lattice, compared with square lattice. 

Now we consider the effect on AF superexchange interaction on spin-spin correla- 
tion. A spin-flux phase lfT0l[16| occurs in the square lattice, where four localized spins 
within each square lie (anti)clockwise at the same plane. In the triangular lattice, as } A f 
increases from to 0.4 (in fact J AF - does not exceed 0.1 in material), the system is in FM 
phase at first and then turns into a specific flux phase and eventually evolves into 120° 
spin configuration phase. For the specific flux phase corresponding to the mediate AF 
superexchange interaction, (S,X ) • S; c with the sites i, j and k belonging to one triangle 
statistically equals to 0.7, and the spin-spin correlation of adjacent sites is —0.3. So this 
phase is different from both FM phase and 120° configuration phase. The peak S(q) is 
located at vectors (^,0), (0,^) and (\,\), respectively, and S(q)/N ~ 0.29 very close to its 
maximum value 1/3 (Fig. |5je)). We further extend the system size from 6x 6 to 12 x 12, 
and find that the peaks position does not change and the normalized peak value S(q)/N 
converges to 0.285 (Fig. [6]). So finite-size effect is small and this flux phase is very stable. 

We systematically investigate the behavior of the spin-spin correlation and spin struc- 
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ture factor at high temperature. The peak value of spin structure factor decreases with the 
temperature, and eventually at beta=10 the system is paramagnetic as Jaf ranges from 
to 0.1. So there is no phase transition at high temperature and the system is still subjected 
to the dimensionality of system. 




Figure 5: The spin structure factor S(qi,cj2) of 6x6 triangular lattice. /3 = 75, and /h = 8. (a) } AF = 0, }i = ~ 6 
and <«>= 0.9272, (b) J AF = Q, fi = -8 and <«> = 0.4512, (c) J AF = 0, fi = -10 and <n>= 0.1944, (d) 
J AF = 0.1, ]i = -b and <n> = 0.9146, (e) / AF = 0.1, pi = -8 and <w> = 0.5, and (f) J AF = 0.l, ^ = -10 and 
<w> = 0.2101. 
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Figure 6: The peak S(q) vs the size L in two dimensional triangular lattices. /3 = 75, /h = 8, J AF = 0.1, }i = —S 
and <n>=0.5. 



4 Conclusion 

Double exchange model with antiferromagnetic spin-spin superexchange interaction has 
been studied for two dimensional triangular lattice with the truncated polynomial ex- 
pansion method. For the half-filling case we obtained 120° spin configuration at low 
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temperature, and it is further stabilized by superexchange ineraction. For the quarter- 
filling case, we found that the FM phase is quite stable, and superexchange interaction 
results in a unique spin-flux configuration with very small finite-size effect. 
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